home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 September / macformat-004.iso / Shareware City / Graphics / VideoToolbox ƒ / VideoToolboxSources / Luminance.c < prev    next >
Encoding:
Text File  |  1994-07-07  |  38.7 KB  |  1,059 lines  |  [TEXT/KAHL]
  1. /*
  2. Luminance.c
  3. See Luminance.h for prototypes.
  4. See Luminance.note for explanation.
  5. Also see:
  6. D.G. Pelli and L. Zhang (1991) Accurate control of contrast on microcomputer displays. 
  7. Vision Research, 31:1337-1360.
  8.  
  9. Copyright (c) 1989-1993 Denis G. Pelli. This software is free; you may use it in
  10. your research and give it away to others, with the following restrictions. Any
  11. copy you give away must include this paragraph, unmodified, and, if you have
  12. made any changes to the rest of the file, must include a note, added to HISTORY
  13. below, giving your name, the date, and a description of the changes. This
  14. software may not be sold, whether in source or compiled form, without my
  15. permission. I hope you will find this software useful, but I can't promise that
  16. it will work for you, and am not offering any support. That's why it's free. I
  17. would appreciate reports of bugs and improvements.
  18.  
  19. Denis G. Pelli, Ph.D.
  20. Institute for Sensory Research, Syracuse University, Syracuse, NY 13244-5290, U.S.A.
  21. denis_pelli@isr.syr.edu
  22.  
  23. HISTORY:
  24.  
  25. 4/24/89 dgp first working version.
  26. 4/28/89 dgp added device argument.
  27. 6/23/89 dgp added support for a video attenuator that has unequal gain in the three
  28.             channels. Added LoadClut() and LuminanceClut() routines.
  29. 7/30/89 dgp replaced call to SetEntries (which calls the Color Manager) by a call
  30.             to GDSetEntries (which calls the device driver).
  31. 8/10/89 dgp fixed bugs.
  32. 8/13/89 dgp rewrote LToV() to use nth order polynomial. Broke up LinearizeClut into
  33.             three routines--SetLuminanceRange, SetLuminance, and SetLuminances--the
  34.             last of which is equivalent to the old LinearizeClut.
  35. 8/16/89    dgp    After reading the Brooktree DAC manual I changed the documentation and
  36.             increased tolerance from 0.5 to 1 LSB, and now use the highest
  37.             luminance gain (worst case) over the requested range to transform it
  38.             into a luminance difference tolerance.
  39. 8/21/89    dgp    Made minor improvements in LToV() to deal with failure to converge due to 
  40.             getting stuck in a local minimum of the quartic. My solution is to now
  41.             require that Lmin and Lmax be filled in by the user in the calibration
  42.             structure. In practice this will guarantee that the minimum luminance
  43.             will not be below the local minimum found at the base of the rising curve.
  44.             Last week I also changed the LToV algorithm slightly, to do a bisection if
  45.             Newton's method would take us out of range.
  46. 9/7/89 dgp    Fixed bug in sorting routine, Sort().
  47. 9/10/89 dgp    Commented out the polynomial versions of VToL and LToV and replace them by
  48.             new versions that use a power law. The power law is a marginally better fit
  49.             than a 6th order polynomial (using only 4 instead of 7 parameters) and its
  50.             inverse can be computed much more quickly, about 0.3 ms instead of 2 ms.
  51. 9/11/89 dgp    Deciding to be indecisive, I introduced a conditional POWER_LAW_FIT
  52.             to control the choice of power law or polynomial fit for LToV and VToL.
  53. 9/16/89    dgp    Having finished writing the paper (Pelli & Zhang, 1991) I came back to
  54.             change the program to agree with what I wrote. I changed the equivalent
  55.             number tolerance to be the sum instead of the maximum of the 
  56.             variable-dac gains.
  57. 9/25/89 dgp Fixed minor bug in computation of tolerance. Now gives correct answer even
  58.             when the lowLuminance or highLuminance is out of range.
  59. 10/9/89 dgp    Made minor change so that when some of the gain[] are zero SetLuminance will
  60.             load zero into the unused colorSpec array elements. This is only a cosmetic
  61.             change.
  62. 10/28/89 dgp Made minor changes. I declared LToVPolynomial and LToVPower and
  63.             VToLPolynomial and VToLPower instead of having a conditional compilation. 
  64.             This makes both flavors of routine always available. I also changed 
  65.             LToVPolynomial to use LToVPower instead of LToVQuadratic for its initial 
  66.             guess. The quadratic guess was often poor and occasionally led to 
  67.             convergence problems. If the power law fit is not available then it reverts
  68.             to using the quadratic fit. If that's not available then it reverts to using
  69.             a middle-of-range guess.
  70. 10/28/89 dgp I eliminated the compile-time flag POWER_LAW_FIT and instead use whichever, 
  71.             polynomial or power, provides a better fit, as determined at run time. I 
  72.             biased the comparison of rms error to favor the power law fit since it has 
  73.             few parameters and can be inverted much more quickly. Note: if the 
  74.             power, polynomial, or quadratic fit parameters are not supplied then the 
  75.             appropriate LR.powerError, LR.polynomialError, or LR.quadraticError field 
  76.             should be set to infinity: INF or 1.0/0.0.
  77. 12/5/89    dgp    Added the trivial routine LToL() which enforces the bounds LP->LMin and 
  78.             LP->LMax.
  79. 4/2/90    dgp    Capitalized the word "To" in all function names, e.g. LToL().
  80. 4/22/90    dgp    Version 1.4. Made minor changes to the documentation. Renamed LToV() to
  81.             LToVFormulaic(), and introduced a new LToV() that uses a table to run faster.
  82.             LToV() takes an average of 170 µs, whereas LToVFormulaic() takes
  83.             1700 µs. I also now use the luminanceTable, if available, to speed up
  84.             VToL() slightly, from 310 µs to 180 µs. Each call to SetLuminance() now
  85.             takes 0.9 ms whereas it used to take 2.5 ms. The loss in accuracy is
  86.             negligible. The interpolation error can be reduced still further by
  87.             increasing LUMINANCES_IN_TABLE. Each doubling of the table size quarters the
  88.             maximum possible error of the interpolation.
  89. 4/23/90    dgp    Added a few lines of code to LToV() to check near the last index, in the
  90.             spirit of Numerical Recipes in C, hunt.c, before starting the bisection
  91.             search. I updated the times in the 4/22/90 note to reflect the latest timing
  92.             by TestLuminance.
  93. 7/27/90    dgp    Tightened up the error checking for illegal entries into the clut. I have
  94.             the impression that, contrary to specification, the Apple video driver
  95.             crashes and corrupts itself if given an out-of-range entry value in a
  96.             setEntry Control call.
  97. 9/18/90    dgp    Changed all instances of "v" to "V". The final version of the Pelli
  98.             & Zhang (1991) paper refers to the nominal voltage v; this file
  99.             now refers to the "equivalent number" V; they are related by V=255*v.
  100. 9/24/90    dgp    Updated the documentation to correspond more closely to the (hopefully)
  101.             final version of the Pelli & Zhang (1991) paper.
  102. 10/29/90 dgp Doubled speed of SetLuminance(). Now SetLuminance() takes 133 ms to
  103.             do a complete clut for a small new luminance range, and takes 188 ms for a
  104.             large new range. This required minor changes to SetLuminanceRange() & LToV().
  105.             Changed all function headers to Standard C prototype style.
  106. 10/30/90 dgp Replaced phrase "clut value" by "nominal voltage", for consistency with
  107.             Pelli & Zhang (1991).
  108.             Introduced new fixed fraction data type called "Milli", defined in Milli.h.
  109.             By judicious replacement of double by Milli, particularly in the luminance
  110.             table, I have increased the speed of SetLuminance by a factor of three.
  111.             A Mac II now takes only 42 to 50 ms to build a whole clut. This new feature
  112.             is enabled by setting FAST_LUMINANCE to 1 in Luminance.h. There is a slight
  113.             increase in the luminance tolerance.
  114. 10/31/90 dgp More fine tuning. Now takes 31 to 44 ms to build a whole clut. Introduced
  115.             conditional FAST_LUMINANCE in Luminance.h that causes the same code to
  116.             be compiled with either Milli or double variables. The _SetLuminance()
  117.             routine is now quite polished, and runs fast. All variables and routines
  118.             whose names begin with underscore _ are written here with type Milli, but
  119.             if FAST_LUMINANCE is false then this type is redefined (in this file only)
  120.             as double, and all the Milli macros are redefined to be appropriate for
  121.             a double. So bear in mind that things beginning with underscore are of
  122.             unknown type. The virtue of being able to run the same code with Milli or
  123.             double computations is that if you suspect an error in the (homemade)
  124.             Milli arithmetic, you can try out double arithmetic simply by setting
  125.             FAST_LUMINANCE to 0 and recompiling.
  126.             In the interest of speed I have streamlined the tolerance
  127.             computation. The new tolerances are probably as good as the old ones.
  128. 11/2/90 dgp Added _SetLuminances() subroutine that substitutes linear interpolation
  129.             for most of the calls to _LToV(). This has greatly speeded up
  130.             SetLuminancesAndRange(), with minor loss of accuracy. A complete clut now
  131.             takes 9 to 31 ms. Even a Mac II can now compute low-contrast cluts
  132.             on the fly, between frames. (The user may wish to optimize the speed-
  133.             accuracy trade-off controlled by the LINEAR_V_DOMAIN parameter, which
  134.             determines the maximum width of the interpolation interval.)
  135. 11/6/90 dgp Replaced Milli by FIXED, which can be compiled as either double or Fixed.
  136.             The Fixed math is sufficiently precise as to give the same DAC values as
  137.             double math.
  138.             Now figure out optimal bit shift LT->L.LShift to preserve as much resolution
  139.             as possible when interpolating in _LToV(). 
  140.             A complete clut now takes 8 to 30 ms. For threshold stimuli the average
  141.             time will usually be less than 15 ms.
  142. 11/8/90 dgp    Eliminated gamma slope table since the speed-up it offered was too small
  143.             to measure. Tidied up the computation of LShift.
  144. 8/24/91    dgp    Made compatible with THINK C 5.0.
  145. 12/17/91 dgp Added new routines: IncrementLuminance() and GetV().
  146.             The former is used to obtain the lowest possible contrast. On a log contrast
  147.             scale, minus infinity (zero contrast) is a poor approximation to even
  148.             a small log contrast.
  149. 3/11/92    dgp    Minor editing of comments.
  150. 12/2/92 dgp Fixed minor bug in SetLuminanceRange() reported by Wei Xie and Ken Alexander
  151.             that could cause the THINK C Debugger to report an error when an 
  152.             uninitialized double V[] was converted to an int. The value was not used 
  153.             for anything, so this was innocuous.
  154. 12/17/92 dgp 1.94 Began enhancement to allow any dacSize, not just 8-bit dacs. I have
  155.             confirmed that everything still works fine with 8-bit dacs, but I don't
  156.             have a 9-bit dac video card to test it any further. •CAUTION: this probably
  157.             does not yet work correctly for 9-bit dacs. E.g. it seems likely that the 
  158.             fixed point math will break if VMax is increased to 511; a lot of careful 
  159.             error analysis went into the original coding, at which time I thought that 
  160.             it was safe to assume an 8-bit dac (i.e. VMax==255). •Removed obsolete 
  161.             support for THINK C 4. 
  162. 12/21/92 dgp No longer load unused dac bits.
  163. 1/4/93    dgp    LoadLuminances now checks the gdType.
  164. 5/10/93    dgp    If driver is in gray mode (i.e. not color) then LoadLuminances maps rgb to 
  165.             gray by simply copying the green component to the red and blue components.
  166. 5/12/93    dgp removed debugging code that forced isGray to always be true.
  167. 3/24/94    dgp    added LToE, LToEOrdered, and EToL.
  168. */
  169. #include "VideoToolbox.h"        /* prototype for GDSetEntries() */
  170. #include "Luminance.h"
  171. #include <assert.h>
  172. #include <math.h>
  173.  
  174. static void Sort(int n,double arr[],int krr[]);    /* for internal use only */
  175.  
  176. #undef LongToFix
  177. #undef FixToLong
  178. #undef DoubleToFix
  179. #undef FixToDouble
  180. #if FAST_LUMINANCE
  181.     #define FIXED Fixed
  182.     #if !MACINTOSH
  183.         typedef long Fixed;
  184.         #define FixMul(x,y) DoubleToFix(FixToDouble(x)*FixToDouble(y))
  185.         #define FixDiv(x,y) DoubleToFix(FixToDouble(x)/FixToDouble(y))
  186.     #endif
  187.     #define LongToFix(x) ((long)(x)<<16)
  188.     #define FixToLong(x) ((x)>>16)
  189.     #define DoubleToFix(x) ((Fixed)((x)*65536.+0.5))
  190.     #define FixToDouble(x) ((double)(x)*(1./65536.))
  191.     #define D(x) FixToDouble(x)                            /* handy for debugging */
  192. #else
  193.     #define FIXED double
  194.     #define FixMul(x,y) ((double)(x)*(y))
  195.     #define FixDiv(x,y) ((double)(x)/(y))
  196.     #define LongToFix(x) ((double)(x))
  197.     #define FixToLong(x) ((long)(x))
  198.     #define DoubleToFix(x) ((double)(x))
  199.     #define FixToDouble(x) ((double)(x))
  200. #endif
  201.  
  202. #undef MAX
  203. #undef MIN
  204. #define MAX(a,b) ((a) > (b) ? (a) : (b))
  205. #define MIN(a,b) ((a) < (b) ? (a) : (b))
  206.  
  207. double EToL(luminanceRecord *LP,int entry)
  208. {
  209.     return GetLuminance(NULL,LP,entry);
  210. }
  211.  
  212. int LToE(luminanceRecord *LP,double L,int firstEntry,int lastEntry)
  213. //Returns the index of the table entry in specified range with luminance closest to L.
  214. {
  215.     int i,entry;
  216.     double dL,latestDL;
  217.     
  218.     latestDL=1.0/0.0;
  219.     entry=0;
  220.     for(i=firstEntry;i<=lastEntry;i++){
  221.         dL=fabs(L-GetLuminance(NULL,LP,i));
  222.         if(dL<latestDL){
  223.             entry=i;
  224.             latestDL=dL;
  225.         }
  226.     }
  227.     return entry;
  228. }
  229.  
  230. int LtoEOrdered(luminanceRecord *LP,double L,int firstEntry,int lastEntry)
  231. /*
  232. Returns the index of the table entry in specified range with luminance closest to L.
  233. Assumes ordered table from firstEntry to lastEntry, either increasing or decreasing.
  234. Returned entry is guaranteed to be in the range firstEntry...lastEntry.
  235. Based on Numerical Recipes "locate.c", page 117
  236. */
  237. {
  238.     register int lower,middle,upper;
  239.     register double LLower,LMiddle,LUpper;
  240.     int ascend;
  241.     
  242.     lower=firstEntry;
  243.     upper=lastEntry;
  244.     LLower=GetLuminance(NULL,LP,lower);
  245.     LUpper=GetLuminance(NULL,LP,upper);
  246.     ascend=(LUpper>LLower);
  247.     while(upper-lower>1){
  248.         middle=(upper+lower)/2;
  249.         LMiddle=GetLuminance(NULL,LP,middle);
  250.         if(L>LMiddle == ascend){
  251.             lower=middle;
  252.             LLower=LMiddle;
  253.         }
  254.         else{
  255.             upper=middle;
  256.             LUpper=LMiddle;
  257.         }
  258.     }
  259.     if(L-LLower<LUpper-L == ascend)return lower;
  260.     else return upper;
  261. }
  262. double SetLuminance(GDHandle device,luminanceRecord *LP
  263.     ,int theEntry,double luminance
  264.     ,double lowLuminance,double highLuminance)
  265. /*
  266. Set one entry in the ColorSpec table (and the clut if device is not NULL) to
  267. a specified luminance. It's ok for lowLuminance to be greater than highLuminance; 
  268. they still designate a range.
  269. */
  270. {
  271.     FIXED _luminance;
  272.     
  273.     if(theEntry<0 || theEntry>=COLORS 
  274.         || device!=NULL && theEntry>=GDCLUTSIZE(device))
  275.             PrintfExit("\007SetLuminance: illegal entry %d\n",theEntry);
  276.     
  277.     SetLuminanceRange(LP,lowLuminance,highLuminance);
  278.     _luminance=DoubleToFix(luminance);
  279.     _SetLuminance(LP,theEntry,_luminance);
  280.     if(device != NULL)LoadLuminances(device,LP,theEntry,theEntry);
  281.     return FixToDouble(_Tolerance(LP,_luminance));
  282. }
  283.  
  284. double SetLuminances(GDHandle device,luminanceRecord *LP
  285.     ,int firstEntry,int lastEntry
  286.     ,double firstLuminance,double lastLuminance)
  287. /*
  288. Set a series of entries in the ColorSpec table (and the clut if device is not NULL)
  289. to a linear sequence of luminances. Assume this is the entire luminance range of 
  290. interest.
  291. */
  292. {
  293.     return SetLuminancesAndRange(device,LP,firstEntry,lastEntry
  294.         ,firstLuminance,lastLuminance,firstLuminance,lastLuminance);
  295. }
  296.  
  297. double SetLuminancesAndRange(GDHandle device,luminanceRecord *LP
  298.     ,int firstEntry,int lastEntry
  299.     ,double firstLuminance,double lastLuminance
  300.     ,double lowLuminance,double highLuminance)
  301. /*
  302. Set a series of entries in the ColorSpec table (and the clut if device is not
  303. NULL) to a linear sequence of luminances. Uses last two arguments to set the
  304. luminance range.
  305.  
  306. I introduced a stack-space check before calling _SetLuminances(), since it calls
  307. itself recursively and may use 1000 bytes all together. However, printing an
  308. error message requires at least 4500 bytes. So I quit if the stack gets below
  309. 6000, so that the user will be presented with an understandable error message
  310. rather than a mysterious quit or crash.
  311. */
  312. {
  313.     FIXED _tolerance,_t,_dL64;
  314.     FIXED _firstL,_lastL,_firstV,_lastV;
  315.     
  316.     if(firstEntry<0 || firstEntry>lastEntry || lastEntry>=COLORS 
  317.         || device!=NULL && lastEntry>=GDCLUTSIZE(device) )
  318.         PrintfExit("\007SetLuminancesAndRange: illegal entries %d %d\n",firstEntry,lastEntry);
  319.     if(StackSpace()<6000)PrintfExit("SetLuminancesAndRange: not enough stack space.\007\n");
  320.     SetLuminanceRange(LP,lowLuminance,highLuminance);
  321.     _firstL=DoubleToFix(firstLuminance);
  322.     _lastL=DoubleToFix(lastLuminance);
  323.     _firstV=_LToV(LP,_firstL + LP->_LOffset);
  324.     _lastV=_LToV(LP,_lastL + LP->_LOffset);
  325.     if(lastEntry != firstEntry){
  326.         #if FAST_LUMINANCE
  327.         _dL64=DoubleToFix(64.0*(lastLuminance-firstLuminance)/(lastEntry-firstEntry));
  328.         #else
  329.         _dL64=DoubleToFix((lastLuminance-firstLuminance)/(lastEntry-firstEntry));
  330.         #endif
  331.     }else _dL64=0;
  332.     _SetLuminances(LP,firstEntry,lastEntry,_firstL,_dL64,_firstV,_lastV);
  333.     if(device != NULL)LoadLuminances(device,LP,firstEntry,lastEntry);    
  334.     
  335.     /* quick and dirty estimate of tolerance */
  336.     _tolerance=LP->L._L[0] - _firstL;
  337.     _t=_lastL - LP->L._L[LUMINANCES_IN_TABLE-1];
  338.     if(_tolerance<_t)_tolerance=_t;
  339.     if(_tolerance<0)_tolerance=0;
  340.     _tolerance+=LP->_tolerance;
  341.     return FixToDouble(_tolerance);    /* estimate of max error in ΔL */
  342. }
  343.  
  344. void _SetLuminances(luminanceRecord *LPtr,int first,int last
  345.     ,FIXED _firstL,FIXED _dL64,FIXED _firstV,FIXED _lastV)
  346. /*
  347. This routine eliminates most of the the slow _LToV() table lookups and
  348. interpolations by assuming that the gamma function is smooth enough that we may
  349. linearly interpolate over any interval of V no wider than LINEAR_V_DOMAIN. Since
  350. the gamma function is roughly parabolic, the consequent error in L will be
  351. proportional to the square of LINEAR_V_DOMAIN, and independent of where this
  352. interval is along V. Since we're inscribing straight line segments in a
  353. positively accelerating gamma function, we will overestimate luminance, and
  354. consequently V will be too low. The error will be greatest at the middle of each
  355. linearly interpolated V interval.
  356.  
  357. If the requested V interval is larger than LINEAR_V_DOMAIN, then it is split
  358. into two intervals by making a pair of recursive calls. LINEAR_V_DOMAIN is
  359. defined in Luminance.h. I suggest a value of 4, but it may be set larger to
  360. attain slightly higher speed at lower accuracy.
  361.  
  362. CAUTION: This routine is a stack hog. Each call uses up about 64 bytes on the
  363. stack, and it calls itself recursively, up to log2(256/LINEAR_V_DOMAIN) times.
  364. Do NOT change any of the declarations of the "new" variables used in the first
  365. if{} to "static", as that would cause the recursion to fail.
  366. */
  367. {
  368.     static int doLast=1;
  369.     int saveDoLast;
  370.     int newOne;
  371.     FIXED _newL,_newV;
  372.     register luminanceRecord *LP=LPtr;
  373.     register int i,Vi;
  374.     register FIXED _VToGo,_g;
  375.     static RGBColor *RGBPtr;    /* static in order to minimize stack usage */
  376.     static int theEntry;        /* static in order to minimize stack usage */
  377.     static FIXED _V,_dV;        /* static in order to minimize stack usage */
  378.     register short leftShift;
  379.     
  380.     if(last-first>1 ){
  381.         _dV=_lastV-_firstV;
  382.         if(_dV>LongToFix(LINEAR_V_DOMAIN) || _dV<LongToFix(-LINEAR_V_DOMAIN)){
  383.             newOne=(first+last)>>1;
  384.             #if FAST_LUMINANCE
  385.                 _newL=_firstL+((newOne-first)*_dL64>>6);
  386.             #else
  387.                 _newL=_firstL+(newOne-first)*_dL64;
  388.             #endif
  389.             _newV=_LToV(LP,_newL + LP->_LOffset);
  390.             saveDoLast=doLast;
  391.             doLast=0;
  392.             _SetLuminances(LP,first,newOne,_firstL,_dL64,_firstV,_newV);
  393.             doLast=saveDoLast;
  394.             _SetLuminances(LP,newOne,last,_newL,_dL64,_newV,_lastV);
  395.             return;
  396.         }
  397.     }
  398.     if(last!=first)_dV=(_lastV-_firstV)/(last-first);
  399.     else _dV=0;
  400.     if(!doLast)last--;
  401.     leftShift=LP->leftShift;
  402.     for(theEntry=first,_V=_firstV;theEntry<=last;theEntry++,_V+=_dV){
  403.         /****** This section of code is copied from _SetLuminance() ****/
  404.         RGBPtr = &LP->table[theEntry].rgb;
  405.         *RGBPtr=LP->rgb;                        /* load for fixed DACs */
  406.         _VToGo=_V - LP->_VFixed + LP->_VHalfStep;
  407.         for(i=LP->fixed;i<LP->dacs;i++) {
  408.             _g=LP->_gain[i];
  409.             Vi=_VToGo/_g;                        /* truncate to integer */
  410.             if(Vi>LP->VMax)Vi=LP->VMax;
  411.             if(Vi<LP->VMin)Vi=LP->VMin;
  412.             _VToGo -= Vi*_g;
  413.             ((short *)RGBPtr)[LP->dac[i]]=Vi<<leftShift;
  414.         }
  415.         /***************************************************************/
  416.     }
  417.     return;
  418. }
  419.  
  420. void _SetLuminance(luminanceRecord *LPtr,int theEntry,FIXED _luminance)
  421. /*
  422. Set one entry in the ColorSpec table to a specified luminance. This is the
  423. private subroutine that actually does all the work for SetLuminance(). This
  424. routine is designed to run as fast as possible.
  425. */
  426. {
  427.     register luminanceRecord *LP=LPtr;
  428.     register int i,Vi;
  429.     register FIXED _VToGo,_g;
  430.     RGBColor *RGBPtr;
  431.     register short leftShift=LP->leftShift;
  432.     
  433.     RGBPtr = &LP->table[theEntry].rgb;
  434.     *RGBPtr=LP->rgb;                        /* load for fixed DACs */
  435.     _VToGo=_LToV(LP,_luminance + LP->_LOffset) - LP->_VFixed + LP->_VHalfStep;
  436.     for(i=LP->fixed;i<LP->dacs;i++) {
  437.         _g=LP->_gain[i];
  438.         Vi=_VToGo/_g;                        /* truncate to integer */
  439.         if(Vi>LP->VMax)Vi=LP->VMax;
  440.         if(Vi<LP->VMin)Vi=LP->VMin;
  441.         _VToGo -= Vi*_g;
  442.         ((short *)RGBPtr)[LP->dac[i]]=Vi<<leftShift;
  443.     }
  444. }
  445.  
  446. void IncrementLuminance(GDHandle device,luminanceRecord *LPtr,int theEntry)
  447. /*
  448. Make smallest possible increase of the luminance of one entry in the ColorSpec
  449. table.
  450. */
  451. {
  452.     double V;
  453.     register luminanceRecord *LP=LPtr;
  454.     
  455.     V=GetV(device,LP,theEntry);
  456.     V+=LP->VHalfStep;
  457.     V+=LP->VHalfStep;
  458.     SetLuminance(device,LP,theEntry,VToL(LP,V),LP->lowLuminance,LP->highLuminance);
  459. }
  460.  
  461. FIXED _Tolerance(luminanceRecord *LPtr,FIXED _luminance)
  462. {
  463.     register luminanceRecord *LP=LPtr;
  464.     register FIXED _tolerance,_t;
  465.     
  466.     /* quick and dirty estimate of tolerance */
  467.     if(LP->L.exists==luminanceSet){
  468.         _tolerance=LP->L._L[0] - _luminance;
  469.         _t=_luminance - LP->L._L[LUMINANCES_IN_TABLE-1];
  470.         if(_tolerance<_t)_tolerance=_t;
  471.         if(_tolerance<0)_tolerance=0;
  472.         _tolerance+=LP->_tolerance;
  473.     }
  474.     else _tolerance=LP->_tolerance;
  475.     return _tolerance;
  476. }
  477.  
  478. double GetLuminance(GDHandle device,luminanceRecord *LP,int theEntry)
  479. /*
  480. If device is not NULL then examines one entry in the actual clut, otherwise
  481. examines the ColorSpec table contained in *LP, and in either case returns the
  482. luminance that will be produced. Supplying an illegal entry value results in a
  483. returned value of -INF.
  484. */
  485. {
  486.     return VToL(LP,GetV(device,LP,theEntry));
  487. }
  488.  
  489. double GetV(GDHandle device,luminanceRecord *LP,int theEntry)
  490. /*
  491. If device is not NULL then examines one entry in the actual clut, otherwise
  492. examines the ColorSpec table contained in *LP, and in either case returns the V
  493. that will be produced. Supplying an illegal entry value results in a returned
  494. value of -INF.
  495. */
  496. {
  497.     RGBColor *myRGBPtr=NULL;
  498.     ColorSpec myCSpec;
  499.     double V;
  500.     int error;
  501.     
  502.     if(device != NULL) {
  503.         if(theEntry<0 || theEntry>=GDCLUTSIZE(device)){
  504.             PrintfExit("GetLuminance: illegal entry %d\n",theEntry);
  505.         }
  506.         error=GDGetEntries(device,theEntry,0,&myCSpec);
  507.         if(error){
  508.             PrintfExit("GetLuminance: GDGetEntries error %d\007",error);
  509.         }
  510.         myRGBPtr=&myCSpec.rgb;
  511.     }
  512.     else {
  513.         if(LP!=NULL && theEntry>=0 && theEntry<COLORS) 
  514.             myRGBPtr=&LP->table[theEntry].rgb;
  515.     }
  516.     if(myRGBPtr == NULL) return -INF;
  517.     V = LP->r*(myRGBPtr->red>>8);
  518.     V += LP->g*(myRGBPtr->green>>8);
  519.     V += LP->b*(myRGBPtr->blue>>8);
  520.     return V;
  521. }
  522.  
  523. void LoadLuminances(GDHandle device, luminanceRecord *LP,
  524.     int firstEntry, int lastEntry)
  525. /*
  526. This just calls GDSetEntries() or GDDirectSetEntries() to load your ColorSpec
  527. table into the clut of your screen device. It is here simply to provide a
  528. cosmetic match to the call to SetLuminances(), for which loading the clut is
  529. optional. Note: if you prefer, instead of LP you may send just the address of a
  530. ColorSpec table, cast to (luminanceRecord *), since a luminanceRecord begins
  531. with a ColorSpec table.
  532. */
  533. {
  534.     int error;
  535.     short isGray,i;
  536.     unsigned short j;
  537.     RGBColor *RGBPtr;
  538.     ColorSpec table[256],*tablePtr;
  539.  
  540.     isGray=!TestDeviceAttribute(device,gdDevType);
  541.     if(isGray){
  542.         for(i=firstEntry;i<=lastEntry;i++){
  543.             RGBPtr=&table[i].rgb;
  544.             RGBPtr->red=RGBPtr->green=RGBPtr->blue=LP->table[i].rgb.green;
  545.         }
  546.         tablePtr=&table[firstEntry];
  547.     }else tablePtr=&LP->table[firstEntry];
  548.  
  549.     switch((*device)->gdType){
  550.     case fixedType:
  551.         printf("LoadLuminances: this device has a fixed CLUT.\n");
  552.         return;
  553.     case clutType:
  554.         error=GDSetEntries(device,firstEntry,lastEntry-firstEntry,tablePtr);
  555.         break;
  556.     case directType:
  557.         error=GDDirectSetEntries(device,firstEntry,lastEntry-firstEntry,tablePtr);
  558.         break;
  559.     }
  560.     if(error) printf("\007LoadLuminances: error %d\n",error);
  561. }
  562.  
  563. double LToL(luminanceRecord *LP,double L)
  564. {
  565.     if(L > LP->LMax) return LP->LMax;
  566.     if(L < LP->LMin) return LP->LMin;
  567.     return L;
  568. }
  569.  
  570. double VToL(luminanceRecord *LP,double V)
  571. /*
  572. Return the luminance that would result from a nominal voltage.
  573. */
  574. {
  575.     return FixToDouble(_VToL(LP,DoubleToFix(V)));
  576. }
  577.  
  578. FIXED _VToL(luminanceRecord *LP,FIXED _V)
  579. /*
  580. Return the luminance that would result from a nominal voltage.
  581. */
  582. {
  583.     register int i;
  584.     register luminanceTable *LT;
  585.     register FIXED _di;
  586.     double LF,VF;
  587.     
  588.     LT=&LP->L;
  589.     if(LT->exists==luminanceSet){
  590.         _di=FixDiv(_V-LT->_VMin,LT->_dV);
  591.         i=FixToLong(_di);
  592.         _di -= LongToFix(i);
  593.         if(i<0)return LT->_L[0];
  594.         if(i>=LUMINANCES_IN_TABLE-1)return LT->_L[LUMINANCES_IN_TABLE-1];
  595.         return LT->_L[i]+FixMul(_di,LT->_L[i+1]-LT->_L[i]);
  596.     }
  597.     VF=FixToDouble(_V);
  598.     if(LP->powerError < 2.0*LP->polynomialError) LF=VToLPower(LP,VF);
  599.     else LF=VToLPolynomial(LP,VF);
  600.     return DoubleToFix(LF);
  601. }
  602.  
  603. double LToV(luminanceRecord *LP,double L)
  604. {
  605.     return FixToDouble(_LToV(LP,DoubleToFix(L)));
  606. }
  607.  
  608. FIXED _LToV(luminanceRecord *LP,FIXED _L)
  609. /*
  610. New, faster version uses a tabulated version of the gamma function VToL(). This
  611. replaces the old function LToVFormulaic(), formerly called LToV(). A subtlety is
  612. that I tabulate the gamma function _L[]=VToL() rather than its inverse, because
  613. this yields a much lower upper bound on the error of the linearly interpolated
  614. L. This is true because the gamma function is roughly parabolic. Equally spaced
  615. samples of a parabolic function yield equal maximum error of linear
  616. interpolation in all intervals. The cost of tabulating _L[] rather than its
  617. inverse is that we must search the table in order to find the largest _L[] less
  618. than or equal to L. The search time is minimized by using a search procedure
  619. that begins the search at the place found by the last search, since successive
  620. calls to LToV() are likely to request similar luminances.
  621. */
  622. {
  623.     register luminanceTable *LT;
  624.     register int lo,hi,i;
  625.     register FIXED _dL,_dLTemp;
  626.     FIXED _V,_LLo;
  627.     
  628.     LT=&LP->L;
  629.     if(LT->exists != luminanceSet){
  630.         /* get domain of the monotonic section of the gamma function */
  631.         LT->_VMin=DoubleToFix(LToVFormulaic(LP,LP->LMin));
  632.         LT->_VMax=DoubleToFix(LToVFormulaic(LP,LP->LMax));
  633.         LT->_dV=(LT->_VMax-LT->_VMin)/(LUMINANCES_IN_TABLE-1);
  634.         _V=LT->_VMin;
  635.         for(i=0;i<LUMINANCES_IN_TABLE;i++,_V+=LT->_dV) LT->_L[i]=_VToL(LP,_V);
  636.         /* just in case, impose monotonicity, from center on out */
  637.         for(i=LUMINANCES_IN_TABLE/2;i>0;i--)
  638.             if(LT->_L[i-1] > LT->_L[i]) LT->_L[i-1]=LT->_L[i];
  639.         for(i=1+LUMINANCES_IN_TABLE/2;i<LUMINANCES_IN_TABLE;i++)
  640.             if(LT->_L[i-1] > LT->_L[i]) LT->_L[i]=LT->_L[i-1];
  641.         #if FAST_LUMINANCE    /* compute the slope dL/dV */
  642.             /* find shift that maximizes ΔL precision without overflow */
  643.             _dL=0;
  644.             for(i=0;i<LUMINANCES_IN_TABLE-1;i++){
  645.                 _dLTemp=LT->_L[i+1]-LT->_L[i];
  646.                 if(_dL<_dLTemp)_dL=_dLTemp;
  647.             }
  648.             LT->LShift=Log2L(LONG_MAX/_dL);
  649.         #endif                
  650.         LT->exists=luminanceSet;
  651.         LT->latestIndex=-2;    /* invalid, so hunt will start from scratch */
  652.     }
  653.     /* hunt for L in LT->_L[] */
  654.     /* first check at latestIndex, latestIndex±1, ... */
  655.     lo=-1;
  656.     hi=LUMINANCES_IN_TABLE;
  657.     i=LT->latestIndex;
  658.     if(i>lo && i<hi){
  659.         if(_L>LT->_L[i])
  660.             {lo=i;i++;}
  661.         else
  662.             {hi=i;i--;}
  663.     }
  664.     /* simple bisection search, see Numerical Recipes in C, pages 98-99. */
  665.     while(hi-lo>1){
  666.         i=(hi+lo)>>1;
  667.         if(_L>LT->_L[i])lo=i;
  668.         else hi=i;
  669.     }
  670.     LT->latestIndex=lo;
  671.     if(lo<0){
  672.         LT->latestIndex=0;            /* nearest legal index */
  673.         return LT->_VMin;
  674.     }
  675.     if(lo>=LUMINANCES_IN_TABLE-1)return LT->_VMax;
  676.     _LLo=LT->_L[lo];
  677.     _dL=_L-_LLo;
  678.     #if FAST_LUMINANCE
  679.         return LT->_VMin+lo*LT->_dV
  680.             +(_dL<<LT->LShift)/((LT->_L[lo+1]-_LLo<<LT->LShift)/LT->_dV);
  681.     #else
  682.         return LT->_VMin+LT->_dV*(lo+_dL/(LT->_L[lo+1]-_LLo));
  683.     #endif
  684. }
  685.  
  686. double LToVFormulaic(luminanceRecord *LP,double L)
  687. /*
  688. Find a nominal voltage that would give luminance L. Failing that, it returns the
  689. closest possible value. The answer is computed by inverting the formula for the
  690. gamma function.
  691. */
  692. {
  693.     if(LP->powerError < 2.0*LP->polynomialError) return LToVPower(LP,L);
  694.     else return LToVPolynomial(LP,L);
  695. }
  696.  
  697. double VToLPolynomial(luminanceRecord *LP,double V)
  698. /*
  699. Return the luminance that would result from a nominal voltage. Uses a polynomial
  700. equation L = p[0] + p[1]*V^1 +... of any degree. dgp 8/11/89
  701. */
  702. {
  703.     register double L,VV;
  704.     register int i,m;
  705.     
  706.     m=LP->coefficients;
  707.     L=0.0;
  708.     VV=1.0;
  709.     for(i=0;i<m;i++){
  710.         L+=LP->p[i]*VV;
  711.         VV*=V;
  712.     }
  713.     return L;
  714. }
  715.  
  716. double VToLPower(luminanceRecord *LP,double V)
  717. /*
  718. Return the luminance that would result from a nominal voltage.
  719. Uses a rectified power law:
  720.  
  721.     L(V)=power[0]+Rectify(power[1]+power[2]*V)^power[3]
  722.     
  723.     Rectify(x)=x if x≥0
  724.     Rectify(x)=0 if x<0
  725.  
  726. dgp 10/28/89
  727. */
  728. {
  729.     double L;
  730.     
  731.     L=LP->power[1]+LP->power[2]*V;
  732.     if(L>0.0) L=LP->power[0]+pow(L,LP->power[3]);
  733.     else L=LP->power[0];
  734.     return L;
  735. }
  736.  
  737. double LToVPolynomial(luminanceRecord *LP,double L)
  738. /*
  739. Find a nominal voltage that would give luminance L. Failing that, it returns the
  740. closest possible value. Solves a polynomial equation L = p[0] + p[1]*V^1 +... of
  741. any degree. First we find a quick approximate answer by solving the power law
  742. fit, LToVPower(). Then we use Newton's method to home in quickly on the solution
  743. to the polynomial fit. Newton's method is taken from Numerical Recipes in C. For
  744. a polynomial of degree 8 (which is what I recommend) this routine now does about
  745. 3? iterations and takes about 1.1? ms. The error in V is less than TOLERANCE,
  746. i.e. relative to full scale of 256 we get an accuracy of 14.6 bits. You can
  747. increase the TOLERANCE to reduce the number of iterations to make it slightly
  748. faster. A TOLERANCE of 1.0 (for 8-bit accuracy) takes 0.8? ms, not much of a
  749. savings in time, and a large loss in accuracy. Reducing TOLERANCE to 1e-8 yields
  750. an accuracy of nearly 35 bits and takes 1.5? ms, only slightly longer.
  751.  
  752. This routine will be called essentially once per clut entry, so computing a
  753. whole new clut with entries 0..255 will take 256*1.1? ms = 282? ms. Some
  754. experiments use only part of the clut, and will therefore take less time. In
  755. some of the experiments that do use the whole clut it may be desirable to
  756. compute the clut table in advance.
  757. dgp 8/11/89
  758. */
  759. {
  760.     #define TOLERANCE 0.01                    /* desired accuracy of V, see above */
  761.     #define JMAX 20
  762.     double a[MAX_COEFFICIENTS];
  763.     register int i;
  764.     int m,j;
  765.     register double f,df,u,VV,V,dV;
  766.     
  767.     if(L>LP->LMax) return LP->VMax;
  768.     if(L<LP->LMin) return LP->VMin;
  769.  
  770.     m=LP->coefficients;
  771.     if(m>MAX_COEFFICIENTS || m<1)
  772.         PrintfExit("LToVPolynomial: %d coefficients is too many or too few\007\n",m);
  773.     for(i=0;i<m;i++) a[i]=LP->p[i];
  774.     a[0]-=L;
  775.     if(LP->powerError < 0.2*LP->LMax) V=LToVPower(LP,L);    /* very good initial guess */
  776.     else {
  777.         if(LP->quadraticError < 0.2*LP->LMax) V=LToVQuadratic(LP,L);/* fair initial guess */
  778.         else PrintfExit("LToVPolynomial: neither power nor quadratic fit is available\007\n");
  779.     }
  780.     for(j=0;j<JMAX;j++){
  781.         /* evaluate function, i.e. the polynomial, and its derivative */
  782.         f=a[0];
  783.         df=0.0;
  784.         VV=1.0;
  785.         for(i=1;i<m;i++){
  786.             u=a[i]*VV;
  787.             df+=i*u;
  788.             f+=V*u;
  789.             VV*=V;
  790.         }
  791.         dV=f/df;    /* apply Newton's method */
  792.         if(V-dV<LP->VMin) dV=(V-LP->VMin)/2.0;    /* bisect */
  793.         if(V-dV>LP->VMax) dV=(V-LP->VMax)/2.0;    /* bisect */
  794.         V -= dV;
  795.         if(fabs(dV) < TOLERANCE) return V;
  796.     }
  797.     printf("LToVPolynomial(%7.2f): Warning, too many iterations. VToL(%5.1f)=%7.2f\n",L,V,VToL(LP,V));
  798.     return V;
  799.     #undef TOLERANCE
  800.     #undef JMAX
  801. }
  802.  
  803. double LToVPower(luminanceRecord *LP,double L)
  804. /*
  805. Find a nominal voltage that would give luminance L. Failing that, it returns the
  806. closest possible value. Solves the rectified power law:
  807.     L(V)=power[0]+Rectify(power[1]+power[2]*V)^power[3]
  808.     Rectify(x)=x if x≥0
  809.     Rectify(x)=0 if x<0
  810. LToVPower takes about 300 microseconds on a Mac II with 8881 arithmetic. This
  811. routine will be called essentially once per clut entry, so computing a whole new
  812. clut with entries 0..255 will take 256*0.3 ms = 77 ms. Some experiments use only
  813. part of the clut, and will therefore take less time. In experiments that do use
  814. the whole clut it may be desirable to compute the clut table in advance.
  815. dgp 10/28/89
  816. */
  817. {
  818.     double V;
  819.     
  820.     if(L<LP->power[0]) L=LP->power[0];
  821.     V=(pow(L-LP->power[0],1.0/LP->power[3])-LP->power[1])/LP->power[2];
  822.     if(V>LP->VMax) V=LP->VMax;
  823.     return V;
  824. }
  825.  
  826. double LToVQuadratic(luminanceRecord *LP,double L)
  827. /*
  828. This is a quick, approximate version of LToV() that is used by LToVPolynomial
  829. for an initial guess if no power law fit is available. Find nominal voltage that
  830. would give luminance L. Since a quadratic fit to the luminance calibration is
  831. only a fair approximation, this routine serves only to find a quick approximate
  832. answer to the root of a higher order polynomial fit. Solves the quadratic
  833. equation by method recommended in Numerical Recipes in C. My choice of root
  834. assumes that L is an increasing function of V over the operating range of the
  835. display. If concave (i.e. q[2]<0) use the larger root. If convex (i.e. q[2]>0)
  836. use the smaller root. In plain English, if the parabola is right side up then we
  837. take the right hand branch. If the parabola is upside down then we take the left
  838. hand branch. dgp 8/12/89
  839. */
  840. {
  841.     register double a,b,c,q;
  842.     double V;
  843.  
  844.     if(L>LP->LMax) return LP->VMax;
  845.     if(L<LP->LMin) return LP->VMin;
  846.     c=LP->q[0]-L;
  847.     b=LP->q[1];
  848.     a=LP->q[2];
  849.     if(a == 0.0) return -c/b;
  850.     if(b<0.0)
  851.         q=-0.5*(b-sqrt(b*b-4.0*a*c));
  852.     else
  853.         q=-0.5*(b+sqrt(b*b-4.0*a*c));
  854.     if(a*b<0.0)
  855.         V=q/a;
  856.     else
  857.         V=c/q;
  858.     return V;
  859. }
  860.  
  861. double SetLuminanceRange(luminanceRecord *LPtr,double lowLuminance,double highLuminance)
  862. /*
  863. The user will never call this routine directly, as it is called automatically by
  864. SetLuminance and SetLuminances. This routine does the hard work of figuring out
  865. which DACs to fix to optimally represent a given luminance range. The goal is to
  866. fix the coarser DACs and vary the finer DACs to represent the image. The CHANGES
  867. in luminance (i.e. the contrast) will have the accuracy of the coarsest of the
  868. variable DACs. To avoid the overhead of figuring all this out every single time
  869. you call SetLuminance(), the results are stored temporarily in your
  870. luminanceRecord, and are only recomputed if you request a different luminance
  871. range.
  872.  
  873. This function checks and returns immediately if you call it with the same
  874. luminance range that it already has stored in your luminanceRecord.
  875. */
  876. {
  877.     register luminanceRecord *LP=LPtr;
  878.     register int i;
  879.     double dL,dL2,LOffset;
  880.     double VSum,VFixed;
  881.     int dacs,fixed,dac[DACS];
  882.     double gain[DACS],V[DACS];
  883.     double VGoal;
  884.     double lowV,highV;
  885.     double tolerance;
  886.     double VHalfStep;
  887.     double dV;
  888.     unsigned short *RGBArray;
  889.     
  890.     /* -1. If necessary, swap low & high */
  891.     if(lowLuminance > highLuminance){
  892.         dL=lowLuminance;
  893.         lowLuminance=highLuminance;
  894.         highLuminance=dL;
  895.     }
  896.  
  897.     /* 0. Return right away if our work has already been done */
  898.     if(LP->rangeSet == luminanceSet &&
  899.         LP->lowLuminance==lowLuminance &&
  900.         LP->highLuminance==highLuminance)
  901.             return LP->tolerance;
  902.  
  903.     /* 1. sort the DACs by gain, discard any with zero or negative gain.
  904.     Gain is defined as the change in V when a single DAC is increased from 0 to 1.
  905.     Sort the gains so that gain[0]≥gain[1]≥gain[2]...
  906.     Note this order is opposite to that used in Pelli & Zhang (1991). The answers
  907.     are the same, but the notation is different.
  908.     */
  909.     if(LP->rangeSet != luminanceSet){
  910.         /* This only needs to be done once, even if the range changes. */
  911.         gain[0]=LP->r;
  912.         gain[1]=LP->g;
  913.         gain[2]=LP->b;
  914.         #if DACS>3
  915.             #error "Current implementation doesn't support more than 3 DACs."
  916.         #endif
  917.         for(i=0;i<DACS;i++)dac[i]=i;
  918.         Sort(DACS,gain,dac);            /* sort so that gain[0]≥gain[1]≥gain[2] */
  919.         dacs=0;
  920.         for(i=0;i<DACS;i++) if(gain[i]>0.0)dacs++;
  921.         VHalfStep=gain[dacs-1]/2.0;        /* half the smallest step */
  922.         for(i=0;i<COLORS;i++)LP->table[i].value=i;
  923.         LP->dacSize=Log2L(LP->VMax*2-1);    // in bits, rounded up
  924.         if(LP->dacSize!=8)printf(
  925.             "Caution: Your video card seems to have %ld-bit video dacs, and Luminance.c\n"
  926.             "has only been tested with 8-bit dacs.\n",LP->dacSize);
  927.         /* Apple's convention is to provide a 16-bit value to the dac. Luminance.c
  928.         saves time and space by computing values with only as much precision as
  929.         needed for the the video card's dac (i.e. 0 to VMax, or dacSize bits). This
  930.         value must be shifted left by leftShift bits to yield a standard 16-bit value.
  931.         Note that Apple normally recommends scaling the value to fill the entire
  932.         16-bit range, e.g. multiplying an 8-bit value by 0x101, but that doesn't
  933.         seem appropriate here, where we know the dac size, and are more concerned with
  934.         step size that with range.
  935.         */
  936.         LP->leftShift=16-LP->dacSize;
  937.     }
  938.     else {
  939.         for(i=0;i<DACS;i++){
  940.             gain[i]=LP->gain[i];
  941.             dac[i]=LP->dac[i];
  942.         }
  943.         dacs=LP->dacs;
  944.         VHalfStep=LP->VHalfStep;
  945.     }
  946.     for(i=0;i<DACS;i++)V[i]=0.0;    /* initialize fixed dac voltages */
  947.     
  948.     /* 2. Transform lowLuminance and highLuminance to lowV and highV. */
  949.     lowV=LToV(LP,lowLuminance);        /* takes 0.15 ms */
  950.     highV=LToV(LP,highLuminance);    /* takes 0.15 ms */
  951.  
  952.     /* 3. Designate the finest dacs as variable until we have enough to cover
  953.     the range, then designate the rest as fixed. (DACs 0..fixed-1 will be fixed.)
  954.     */
  955.     VGoal=fabs(highV-lowV);
  956.     VSum=0.0;
  957.     fixed=0;
  958.     for(i=dacs-1;i>=0;i--) {
  959.         if(VSum >= VGoal)fixed++;
  960.         VSum += gain[i]*(LP->VMax - LP->VMin);
  961.     }
  962.     /* sum gains of all variable dacs. This is called gVary by Pelli & Zhang (1991) */
  963.     tolerance=0.0;
  964.     for(i=dacs-1;i>=fixed;i--)tolerance+=gain[i];
  965.     /* scale by highest luminance gain in the requested range */
  966.     dL=fabs(VToL(LP,lowV+1.0)-VToL(LP,lowV));
  967.     dL2=fabs(VToL(LP,highV)-VToL(LP,highV-1.0));
  968.     tolerance *= MAX(dL,dL2);
  969.  
  970.     /* 4. Temporarily set the variable DACs to their midpoints. Now set the fixed DACs
  971.     to most accurately represent (lowV+highV)/2.
  972.     */
  973.     VGoal=(lowV+highV)/2.0;
  974.     VSum=0.0;
  975.     VFixed=(LP->VMax + LP->VMin)/2.0;
  976.     for(i=fixed;i<dacs;i++) VSum+=gain[i]*VFixed;    /* The mid-voltage of the variable DACs */
  977.     for(i=0;i<fixed;i++) {
  978.         V[i]=floor((VHalfStep+VGoal-VSum)/gain[i]);    /* the unattenuated voltage of DAC i */
  979.         V[i]=MAX(LP->VMin,MIN(LP->VMax,V[i]));
  980.         VSum += V[i]*gain[i];
  981.     }
  982.     VSum=0.0;
  983.     for(i=0;i<fixed;i++) VSum+=gain[i]*V[i];
  984.     VFixed=VSum;
  985.     
  986.     /* 5. The limited precision of the fixed DACs may result in a 
  987.     small offset between the requested and now available range. The offset will be
  988.     at most half a step of the finest of the fixed DACs. It seems reasonable
  989.     to offset the requested luminance range by up to that small amount to better fit
  990.     the now available range. */
  991.     LOffset=0.0;
  992.     if(fixed>0){
  993.         VSum=VFixed;
  994.         for(i=fixed;i<dacs;i++) VSum += LP->VMin*gain[i];
  995.         dV= VSum - MIN(lowV,highV);
  996.         dV= MIN(dV,gain[fixed-1]/2.0);
  997.         if(dV>0.0) LOffset+=VToL(LP,VSum+dV)-VToL(LP,VSum);
  998.         VSum=VFixed;
  999.         for(i=fixed;i<dacs;i++) VSum += LP->VMax*gain[i];
  1000.         dV= MAX(lowV,highV) - VSum;
  1001.         dV= MIN(dV,gain[fixed-1]/2.0);
  1002.         if(dV>0.0) LOffset-=VToL(LP,VSum+dV)-VToL(LP,VSum);
  1003.     }
  1004.     /* Now save this information in the luminanceRecord for future use */
  1005.     for(i=0;i<DACS;i++){
  1006.         LP->dac[i]=dac[i];
  1007.         LP->gain[i]=gain[i];
  1008.     }
  1009.     LP->rangeSet=luminanceSet;
  1010.     LP->lowLuminance=lowLuminance;
  1011.     LP->highLuminance=highLuminance;
  1012.     LP->dacs=dacs;
  1013.     LP->VHalfStep=VHalfStep;
  1014.     LP->fixed=fixed;
  1015.     LP->VFixed=VFixed;
  1016.     LP->tolerance=tolerance;
  1017.     LP->LOffset=LOffset;
  1018.  
  1019.     /* Copy into FIXED variables */
  1020.     for(i=0;i<DACS;i++)LP->_gain[i]=DoubleToFix(LP->gain[i]);
  1021.     LP->_VHalfStep=DoubleToFix(VHalfStep);
  1022.     LP->_VFixed=DoubleToFix(VFixed);
  1023.     LP->_tolerance=DoubleToFix(tolerance);
  1024.     LP->_LOffset=DoubleToFix(LOffset);
  1025.  
  1026.     /* cache the fixed DAC values, and zero the rest for good measure */
  1027.     RGBArray = (unsigned short *) &LP->rgb;    /* treat as an array */
  1028.     for(i=0;i<DACS;i++){
  1029.         RGBArray[LP->dac[i]] = (short)V[i]<<LP->leftShift;
  1030.     }
  1031.     return tolerance;
  1032. }
  1033.  
  1034. static void Sort(int n,double arr[],int krr[])
  1035. /*
  1036. Slightly modified sort routine piksr2() from Numerical Recipes in C. I changed
  1037. "float" to "double", and I changed second array to int. I reversed the order, so
  1038. largest element would be first. I changed it to use conventional c arrays,
  1039. starting at index 0.
  1040. */
  1041. {
  1042.     register int i,j;
  1043.     double a;
  1044.     int k;
  1045.  
  1046.     for(j=1;j<n;j++) {
  1047.         a=arr[j];
  1048.         k=krr[j];
  1049.         i=j-1;
  1050.         while (i >= 0 && arr[i] < a) {
  1051.             arr[i+1]=arr[i];
  1052.             krr[i+1]=krr[i];
  1053.             i--;
  1054.         }
  1055.         arr[i+1]=a;
  1056.         krr[i+1]=k;
  1057.     }
  1058. }
  1059.